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ABSTRACT 

Molecular studies on donkey mitochondrial 
sequences have clearly defined two distinct maternal 
lineages involved in domestication. However, 


domestication histories of these two lineages remain 
enigmatic. We therefore compared several 
population characteristics between these two 
lineages based on global sampling, which included 
171 sequences obtained in this study (including 
Middle Asian, East Asian, and African samples) plus 
536 published sequences (including European, 
Asian, and African samples). The two lineages were 
clearly separated from each other based on whole 
mitochondrial genomes and partial non-coding 
displacement loop (D-loop) sequences, respectively. 
The Clade | lineage experienced an increase in 
population size more than 8 000 years ago and 
shows a complex haplotype network. In contrast, the 
population size of the Clade Il lineage has remained 
relatively constant, with a simpler haplotype network. 
Although the distribution of the two lineages was 
almost equal across the Eurasian mainland, they still 
presented discernible but complex geographic bias 
in most parts of Africa, which are known as their 
domestication sites. Donkeys from sub-Saharan 
Africa tended to descend from the Clade | lineage, 
whereas the Clade II lineage was dominant along the 
East and North coasts of Africa. Furthermore, the 
migration routes inferred from diversity decay 
suggested different expansion across China between 
the two lineages. Altogether, these differences 
indicated non-simultaneous domestication of the two 
lineages, which was possibly influenced by the 
response of pastoralists to the desertification of the 
Sahara and by the social expansion and trade of 
ancient humans in Northeast Africa, respectively. 
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INTRODUCTION 


Unlike other species with a similar historical function (e.g., 
horses), domestic donkeys are underrepresented in the 
scientific literature (Blench, 2000). Importantly, donkeys 
remain an essential means of transport in modern society for 
people living in mountain areas, deserts, and poorer regions of 
the world (Smith & Pearson, 2005; Starkey, 2000). Both 
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molecular data and archaeological evidence strongly support 
an African origin for the domestic donkey (Equus asinus) 
(Beja-Pereira et al., 2004; Kimura et al., 2011; Rossel et al., 
2008). This domestication, together with that of horses, is 
believed to have contributed to mobile pastoralism as well as 
the establishment of ancient overland trade routes and spread 
of ancient civilizations (Rosenbom et al., 2015). Previous 
studies on mitochondrial DNA (mtDNA) have revealed that 
modern donkeys can be clearly separated into two distinct 
clades: one clustered directly with the Nubian wild ass (E. 
africanus africanus), hereafter termed the Clade | lineage; 
another derived from a probably extinct wild ass close to the 
Somali wild ass (E. africanus somaliensis), hereafter termed 
the Clade II lineage (Beja-Pereira et al., 2004; Chen et al., 
2006; Kimura et al., 2011). Based on archeological evidence, 
the "Egyptian hypothesis" proposes that villagers in the Nile 
Valley domesticated the resident Nubian wild ass 
approximately 5 000-6 000 years ago (Clutton-Brock, 1992; 
Epstein, 1971); in contrast, the "pastoralist hypothesis", 
supported by ethnographic, climatic, and linguistic data, states 
that pastoralists in northeastern Africa 6 500—7 000 years ago 
domesticated donkeys in response to the increasing aridity of 
the Sahara (Marshall, 2007). Additionally, the very high 
genetic diversity found in the Arabian Peninsula, as based on 
nuclear microsatellites, indicates an alternative domestication 
center or melting pot (ancient trade areas or routes) 
(Rosenbom et al., 2015). Ancient DNA also illustrates the 
extensive distribution of the Nubian and Somali wild ass in 
Northeast Africa (Kimura et al., 2011), alluding to potential 
geographic overlaps between the Nubian wild ass and 
ancestors of the Clade Il donkeys. As such, whether these two 
distinct lineages were domesticated simultaneously or not 
remains controversial, even though they probably originated 
from two mitochondrially distinct wild asses (Jordana et al., 
2016; Xia et al., 2019). 

If the two lineages were domesticated simultaneously, it is 
reasonable to assume the existence of similar genetic 
dynamics corresponding to their rapid expansion with innate 
mobile characteristics. Therefore, in the present paper, we 
investigated the genetic characteristics of 84 mtDNA genomes 
and 707 D-loop sequences to compare the dynamic 
signatures between these two clades that correspond to 
migration history. Our results will help to further understand 
mtDNA diversity in African donkeys and also provide novel 
insights into their domestication histories. 


MATERIALS AND METHODS 
Sample information 


Blood samples were collected and sequenced from animals 
living in villages or donkey breeding farms all over the world. 


This study was approved by the Ethics and Experimental 
Animal Committee of the Kunming Institute of Zoology, 
Chinese Academy of Sciences, China. We selected samples 
without any discrimination toward age, sex, or location. More 
detailed information, including geographical location and the 
US National Center for Biotechnology Information (NCBI) 
GenBank accession Nos., can be found in Supplementary 
Table S1. Specifically, 5 mL of blood was drawn from healthy 
animals by a veterinarian (sampling picture in Supplementary 
Figure S1) with approval from the above-mentioned Ethics 
and Experimental Animal Committee. The blood samples were 
preserved in 95% alcohol for further genomic analysis. Whole 
DNA was isolated from blood using the standard phenol- 
chloroform method (Ma, 2018) and subsequently quantified 
with NanoDrop 2000 (Thermo Fisher Scientific, USA). 

We obtained a 371 bp length polymerase chain reaction 
(PCR) product within the D-loop region for 94 donkeys 
(Tadzhikistan: n=19; Kyrgyzstan: n=5; Kenya: n=16; Nigeria: 
n-22; Iran: n=30; China: n=2) using primers designed in 
previous research (Kefena et al., 2014) (5-DONK-F: CCC 
AAGGACTATCAAGGAAG-3';DONK-R:5'-GGAATGGCCCTGA 
AGAAAG-3'). We strictly followed the reaction system as 
described in Kefena et al. (2014). Specifically, PCR was 
carried out in a 10 pL reaction volume containing 1 uL of DNA 
(50 ug/uL), 1 uL 10xbuffer (10xTakara Taq™ Buffer Mg? 
plus), 0.2 mmol/L of each dNTP, 1 umol/L of each primer, and 
0.2 U/tube Taq (Taq DNA Polymerase, Takara) under the 
following conditions: initial denaturation at 94 °C for 15 min 
followed by 45 cycles of denaturation each at 94 °C for 1 min, 
hybridization at 56 °C for 1 min, extension at 72 °C for 1 min, 
and final extension at 72 °C for 20 min. The PCR products 
were gel-purified as described in Ma (2018) and then 
sequenced by single-strand PCR using the ABI PRISM™ Dye 
Terminator Cycle Sequencing kit following protocols 
recommended by the manufacturer. Sequencing was 
implemented using ABI-PRISM 3730 standard conditions 
(Applied Biosystems, 2009). For each PCR product, 
sequences were determined in both forward and reverse 
directions for all nucleotide positions to avoid possible artificial 
variations. 

Similarly, we obtained the mtDNA genomes of 56 Chinese 
donkeys and two wild asses using ABI-PRISM 3730 with 
primers designed in this study (see Supplementary Table S2). 
Other mtDNA genomes of 17 donkeys (Tadzhikistan: n=6; 
Kyrgyzstan: n=5; Kenya: n=4; Nigeria: n=2) and two 
Tadzhikistan horses were yielded through lon Torrent 
technology, as per previous study (Chen et al., 2016). Briefly, 
we amplified PCR fragments, constructed libraries, performed 
sequencing using next-generation technology, and de novo 
assembled the reads. For quality control, we followed the 
caveats mentioned in previous mtDNA genome study of 
domestic animals (Shi et al., 2014). Variants that differentiated 
from the GenBank reference sequence under accession No. 
NC_001788 (Xu et al., 1996) were scored. We then manually 
checked the bam file exported by Torrent Suite 5.0.2 to 
confirm the scored variants using Integrative Genomics 


Viewer (Thorvaldsdóttir et al., 2013). 

Collectively, we obtained 94 sequences with a 371 bp 
length restricted to the D-loop region and 77 almost complete 
mitochondrial genomes. Sequences generated in this study 
were deposited in the NCBI archive (GenBank accession 
Nos.: MK650231—MK650286 for 58 mtDNA genomes obtained 
by 3730 sequencer; MK619357-MK619411 for 55 D-loop 
haplotypes of Equus yield in this study; 
SAMN!11432793-SAMN!11432809 and SRX5702272- 
SRX5702273 for another 19 mtDNA genomes sequenced 
through lon Torrent technology. For specified topology and 
wider coverage in geographic structure analysis, we combined 
an additional seven mitochondrial genomes and 536 D-loop 
sequences downloaded from NCBI website 
(http://www.ncbi.nlm.nih.gov/), with GenBank accession Nos. 
listed in Supplementary Table S1. Detailed feature information 
of all 707 sequences are included in Supplementary Table S1. 


Phylogenetic construction and estimation of genetic 
parameters 

We reconstructed the phylogenetic tree of equid mtDNA 
sequences using neighbor-joining and maximum likelihood 
analyses with MEGA 6 software (Tamura et al., 2013). For 
mtDNA genomes, 1 000 bootstrap replicates were conducted 
with the Kimura 2-parameter model, assuming Gamma 
distributed mutation rates for neighbor-joining (NJ) analysis 
and Gamma distributed with invariant sites (G+l) distributed 
mutation rates for maximum likelihood (ML) analysis. 
Additionally, we constructed a Bayesian tree using Bayesian 
evolutionary analysis by sampling trees (BEAST) v1.6.1 
(Drummond & Rambaut, 2007) for phylogenetic confirmation 
with the Gamma distributed HKY substitution model under a 
strict clock rate, assuming constant size coalescence for tree 
prior. The non-coding region from 16 129 nt to 16 360 nt was 
excluded due to the short tandem repeats in both horse and 
donkey (repeat unit: GTGCACCT in horse; CACACCCACAC 
ACCCATGCGCGCA in donkey). 

We truncated the 77 mitochondrial genomes into the 371 bp 
length of the D-loop region to ensure combined analysis with 
the 94 D-loop sequences, which yielded 55 haplotypes. The 
downloaded 536 D-loop sequences were from different 
research, which covered different fragments within the D-loop 
region. Therefore, it was inevitable that a very short 
overlapping fragment (235 bp) would be yielded as we 
retrieved as many samples as possible. For genetic diversity 
analysis, because the downloaded 536 D-loop sequences 
were haplotype data, we selected only haplotypes within each 
area from our sequenced samples for parallel comparison. A 
reduced median-joining network was generated using 
NETWORK  (Bandelt et al, 1999). Average pairwise 
differences were estimated using Arlequin v3.5 (Excoffier & 
Lischer, 2010) The nucleotide diversity (rr), mismatch 
distribution, Tajima's D test, and Fu's Fs test were calculated 
using DnaSP v6 (Rozas et al., 2017). 

We assessed the population dynamics for the mtDNA 
genomes of the two lineages using Bayesian Skyline Plot 
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(BSP) implemented in BEAST (Drummond & Rambaut, 2007). 
We estimated the evolution rate using the previously 
estimated divergence time of 2 million years between the 
horse and donkey from paleontological data (Lindsay et al., 
1980). Migration time to regions distant from the domestication 
center was inferred by estimating the evolution time of 
haplotypes only found in the region to the major haplotype. 

A geographic distribution map was plotted using the 
"rworldmap" package (South, 2011) implemented in R (R- 
Core-Team, 2019). 


RESULTS AND DISCUSSION 


Two clearly separated domesticated donkey clades 

A clearly defined phylogenetic relationship can allow 
comparison of genetic characters between two lineages. Thus, 
we constructed a highly supported phylogenetic tree based on 
the 77 mtDNA genome sequences from the current study 
(including 73 donkeys from China, Tadzhikistan, Kyrgyzstan, 
Kenya, and Nigeria, two Asiatic wild asses, and two 
Tadzhikistan indigenous horses) combined with seven mtDNA 
genomes downloaded from NCBI (including two donkeys, two 
Somali wild asses, and two Asian wild asses) (Supplementary 
Table S1). Using the horse as an outgroup, the Somali wild 
ass formed a sister clade with all donkeys, a topology highly 
differentiated from that based on D-loop sequences, which 
was highly supported by the neighbor-joining, maximum 
likelihood, and Bayesian methods analyzed here (Figure 1). 
This topology was also consistent with previous study that 
concentrated on the major mitochondrial coding regions and 
placed the Somali wild ass as a sister clade outside all 
domestic donkeys with high confidence, although they did not 
label the lineage information (Sun et al., 2016). Here, we 
obtained consistent topologies for the major clades according 
to both whole mitochondrial genomes and coding regions only 
(Figure 1; Supplementary Figure S2). With a view that coding 
regions experience restricted selective sweep, whereas 
selective force is potentially relaxed on the D-loop regions 
(Endicott et al., 2009), it is reasonable to assume that the 
phylogenetic relationship revealed by mitochondrial genome 
data would be much closer to reality. 

The two donkey lineages were clearly separated from each 
other, with high support (Figure 1), further confirming the 
existence of two lineages. Considering the lack of geographic 
representation from mtDNA genomes (58 of 75 domestic 
donkeys were from China), we reconstructed the phylogenetic 
tree using both neighbor-joining and maximum likelihood 
methods with D-loop sequences from the 171 sequences 
obtained in this study and 536 sequences publicly available on 
NCBI (see sample information in Supplementary Table S1 and 
phylogenetic tree in Supplementary Figure S3). Similarly, the 
donkeys were divided into two groups. Although the clades 
showed lower support, which is common in analysis of short 
segments like the D-loop region, all sequences with published 
lineage information clustered into their corresponding clades. 
In addition, the 73 donkey D-loop sequences obtained in this 


54 www.zoores.ac.cn 


study showed consistent lineage sorting, as observed in the 
whole mtDNA genome phylogeny, further supporting the D- 
loop region phylogenetic relationship. Our phylogenetic tree 
based on D-loop sequences also revealed a sister relationship 
between the Somali wild ass and all domesticated donkeys 
(Supplementary Figure S3), consistent with the topology 
revealed by the mitochondrial genome, but inconsistent with 
previous studies showing different topology using D-loop 
sequences (Beja-Pereira et al., 2004; Han et al., 2014). The 
number of haplotypes used for phylogenetic construction may 
account for this discrepancy, as more haplotypes represent 
greater genetic diversity, and subsequently more statistical 
credibility. To test this, we randomly selected donkey 
haplotypes of equal sample size, as previously described in 
Han et al. (2014), and constructed a phylogenetic tree with the 
same parameters in this study. We obtained similar 
phylogenetic topology (Supplementary Figure S4), showing a 
closer relationship shared by the Somali wild ass and 
domesticated Clade II lineage. Therefore, a large sample size 
for phylogenetic analysis can compensate for the limited 
information obtained from short segments. 


Different demographic dynamics of two lineages 

As the donkey samples were credibly sorted, we subsequently 
compared several population characteristics between the two 
lineages. First, we scanned sequence variation within the 
mitochondrial genome for each lineage. The Clade | lineage 
showed a total of 288 varied sites: 258 in the coding region 
and 30 in the non-coding regions (D-loop). The Clade Il 
lineage showed a total of 188 varied sites: 168 in the coding 
region and 20 in the non-coding regions (D-loop). Figure 2A 
illustrates the distribution of nucleotide diversity (rr) within 
each lineage along the genome based on assessment of 200 
nt windows (step size-100 nt) centered at the midpoint. 
Although the highest diversity was observed in the D-loop 
region presenting the short tandem repeats of 
CACACCCACACACCCATGCGCGCA (donkey reference 
NC 001788: from 16 173 nt to 16 340 nt) in both lineages, the 
Clade | lineage possessed significantly more segregation sites 
(average segregation sites per window=1.300) in the coding 
region than the Clade II lineage (average segregation sites per 
window=1.091) (Wilcox test: P=0.016), whereas no significant 
differences were found in segregation counts in the non- 
coding region in the two lineages (average segregation sites 
per window in Clade | and Clade Il were 3.057 and 2.303, 
respectively; Wilcox test: P=0.628). This discrepancy may be 
due to the saturation effect in the  D-loop region, 
demographics, selective forces, and/or other factors 
suggestive of potential independent domestic histories 
between the two lineages. To further address this issue, we 
used mtDNA genomes to reconstruct ancestral population 
dynamics with BSP for each lineage (Figure 2B). A similar 
demographic history would be expected under the assumption 
of simultaneous domestication. However, we observed a 
constant effective population size in the Clade ll lineage 
during most history, whereas the Clade | lineage experienced 
an apparent population increase approximately 8 000 years 
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Figure 1 mtDNA genome phylogeny of Equus 
A: Neighbor-joining tree; B: Maximum likelihood tree; C: Bayesian tree. Colored branches denote different categories. Red: Horse; Green: Asiatic 
wild ass; Yellow: Somali wild ass; Blue: Clade | lineage; Magenta: Clade II lineage. 
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ago, coinciding with the archeological date 7 770+95 before | lineage and double peaks in the Clade II lineage (Figure 2C). 
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Harpending, 1992), Fu's Fs test, and Tajima's D test. The Clade Il lineage, confirming the demographic dynamics 
nucleotide mismatch curve showed a single peak in the Clade revealed by mitochondrial genomes. 
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Figure 2 Population characteristics between Clade | and Clade Il lineages 

A: Sequence variation of mtDNA genome within each lineage. Nucleotide diversity (rr) using 200 nt windows (step size=100 nt) centered at 
midpoint. Schematic (linearized) genetic map of mitochondrial genome is presented. B: Demographic histories for each lineage estimated using 
Bayesian skyline plot from mitochondrial genomes. Shadow means 95% highest posterior density intervals of effective number of females. C: 
Mismatch distribution for each mtDNA lineage. D: Median-joining network based on D-loop region. Red diamonds are haplotypes of horse outgroup. 
Green circles are haplotypes of Asiatic wild ass. Dark red circles are haplotypes of Somali wild ass. Indigo circles are haplotypes of Clade | lineage. 
Magenta circles are haplotypes of Clade II lineage. Circle size is proportional to frequencies of haplotypes. 


Discrepancies of haplotype structure showed similar evolution patterns between the two lineages. 
In addition to demographic estimation, we also investigated We first constructed a reduced median network based on D- 
several other genetic characters to assess whether they loop sequences (Figure 2D). Similar to previous study (Beja- 
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Pereira et al., 2004; Kimura et al, 2011), the Clade Il 
haplotypes were mainly derived from a single major haplotype 
(H16), with a simple star-like shape, whereas the genetic 
architecture of the Clade | haplotypes was more complicated, 
with more universally occurring haplotypes (e.g., H20, H52, 
H22). This was consistent with the much higher genetic 
distance and nucleotide diversity within the Clade | lineage 
(average pairwise distance=1.932 7, SD=1.100 7; Pi=0.008 4, 
SD=0.000 3) than within the Clade II lineage (average 
pairwise distance=1.074 4, SD=0.711 2; Pi=0.0046, SD- 
0.000 5), implying that the Clade | lineage involved many more 
individuals at the beginning of domestication compared with 
the Clade II lineage (Table 1). 


Table 1 Genetic distance and nucleotide diversity of two clades 

Pairwise distance SD (distance) Pi SD (Pi 
1.932 7 1.100 7 0.008 4 0.000 3 
1.0744 0.7112 0.004 6 0.000 5 





Clade | lineage 
Clade II lineage 





Biased distribution of two lineages in Africa 

Independent migration events may result in a geographical 
structure. Therefore, we estimated the migration routes of the 
two lineages by referring to their proportion across the world. 
D-loop sequences sampled from the Balkans and 
microsatellites sampled from northeast Africa, the Near East, 
and the Arabian Peninsula indicated the absence of a 
geographical structure (Pérez-Pardal et al., 2014; Rosenbom 
et al., 2015). Indeed, samples collected from Middle Asia 
(Pakistan, Kazakhstan) and major areas of China 
demonstrated an almost equal proportion of the two lineages 
(Figure 3A). The Arabian Peninsula is assumed to be the 
melting pot from where domestic donkeys migrated to the 
world (Rosenbom et al., 2015). As this area showed a nearly 
equal proportion of the two lineages, it is reasonable that a 
similar pattern is commonly observed across the Eurasian 
mainland (Beja-Pereira et al., 2004). Nevertheless, all 20 
samples collected from Nigeria belonged to the Clade | 
lineage. When we focused on the distribution in Africa, an 
apparent spatial structure was detected: donkeys from sub- 
Saharan Africa (e.g., Ghana, Guinea, Benin, Mali, Senegal, 
South Africa, and Burkina Faso) tended to be descended from 
the Clade | lineage, whereas the Clade Il lineage was 
dominant along the East (e.g., Eritrea, Somalia, Swaziland, 
and Zambia) and North coasts (e.g., Libya, Tunisia, and 
Morocco) (Figure 3A). 

The time of migration to areas distant from the 
domestication center can be inferred by dating the time to the 
major haplotype for the derived haplotypes only found in that 
area. In this way, we inferred that the Clade | lineage migrated 
into sub-Sahara 4 874.2841 817.92 years ago, around the 
commencement of desertification in the Sahara (5 000 to 
7 000 BP) (Marshall, 2000). Due to the excellent tolerance of 
donkeys for deserts, it is reasonable to assume that donkeys, 
rather than horses, were the major means of transport across 
the Sahara during the initial period of desertification. Given the 
biased Clade | lineage distribution in the sub-Sahara, this 


migration time provides possible evidence for the "pastoralist 
hypothesis": i. e., pastoralists in northeastern Africa 
domesticated Clade | lineage donkeys in response to the 
increasing aridity in the Sahara. 

Donkeys are thought to be have been brought into Europe 
by the second millennium BC, possibly through viticulture 
introduction, as the donkey is associated with the Syrian god 
of wine, Dionysus (Meutchieye et al., 2017). Interestingly, the 
majority of sequences from the Iberian Peninsula belonged to 
the Clade II lineage, much different from those collected in 
other parts of Europe (Figure 3A). Additionally, the estimated 
time of arrival in the Iberian Peninsula was 5 336.841 652.56 
years ago, much earlier than the known history of ~4 000 
years (Cardoso et al., 2013). Considering the dominance of 
the Clade II lineage in Morocco, it could be assumed that 
donkeys migrated into the Iberian Peninsula directly through 
the Strait of Gibraltar. Furthermore, the unequal distribution of 
Clade | and Clade || in Europe may account for the 
pronounced footprint in American following the complex 
process of colonization, where an apparent geographic 
structure has been observed (Jordana et al., 2016; Xia et al., 
2019). The Clade II lineage also migrated dominantly along 
East Africa, implying a potential association with ancient social 
expansion, such as the Bantu expansion (Hiernaux, 1968; 
Herrera & Garcia-Bertrand, 2018), and ancient trade routes of 
empires such as the Punt and Aksum (Andrews, 2017; Mark, 
2011; Mukhtar, 1981). 


Different routes of two lineages during expansion to 
China 

As populations expand from centers of origin, genetic diversity 
is lost as a consequence of the limited numbers of individuals 
involved in the expansionist movement. Although no apparent 
geographic structure was detected across China, where the 
two lineages show an almost equal proportion (Figure 3A), 
discrepancies in diversity decay may exist between the two 
lineages if non-simultaneous expansion into China occurred, 
thereby suggesting possible private migration routes. The 
Xinjang Province may be a transportation center as it 
demonstrated almost the highest nucleotide diversity for both 
lineages (Figure 3B, C), consistent with previous studies on 
genetic diversity among various Chinese breeds (Ge et al., 
2007) and written records (Xie, 1987). The genetic diversity of 
Clade | remained relatively high in the Qinghai and Henan 
provinces, but declined gradually towards north and south, 
respectively (Figure 3B), consistent with a migration route 
from Xinjiang to the Guanzhong Plain through the Ningxia and 
Gansu provinces, and finally other areas of China inferred 
previously (Ge et al., 2007). The genetic diversity of Clade II 
remained relatively high in the Inner Mongolia and Yunnan 
provinces, but declined substantially in the middle region of 
China (Figure 3C), consistent with previously inferred 
migration routes from Xinjiang to Inner Mongolia (north 
forward) and to Yunnan (south forward) (Ge et al., 2007). 
These potential migration routes partially match written 
records, which suggest that the "Taihang donkey" was 
introduced into Hebei Province through Inner Mongolia; the 
"Yunnan donkey" was introduced from Xinjiang to Yunnan, as 
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Xinjiang 


Qinghai 


Figure 3 Worldwide geographic distribution of two lineages plotted using “rworldmap” package in R (A) and genetic diversity of Clade | 


(B) and Clade II (C) lineages across China 


A: Blue: Clade | lineage; Magenta: Clade II lineage. B: Dark to light indigo indicates diversity decay. Gray regions indicate no samples available. C: 


Dark to light magenta indicates diversity decay. Gray regions indicate no samples available. 


it shares many morphological traits with the Xinjiang donkey; 
and many other breeds may have been formed through 
migration routes from Xinjiang to Gansu and then Guanzhong 
Plain, Henan, and other areas of China (Xie, 1987). Therefore, 
it is likely that the two lineages expanded into China through 
different migration events, although from the same 
transportation center (Xinjiang). 


CONCLUSIONS 


Overall, our results on the demographic dynamics, haplotype 
structure, diversity decay pattern, and distribution bias in 
Africa are in accordance with non-simultaneous domestication 
events and an independent migration history. Our findings 
revealed that domestication of the donkey may have been 
driven by the response of pastoralists to the desertification of 
the Sahara and by social expansion and trade of ancient 
humans. Future research on population genomes will be 
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needed to increase our 
domestication history. 


understanding of donkey 
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